In [1]:
?eps
Out[1]:
In [2]:
?realmax
Out[2]:
In [3]:
subtypes(AbstractFloat)
Out[3]:
In [4]:
# Default values are for Float64
eps(), realmax(), realmin()
Out[4]:
In [5]:
T=Float32
eps(T), realmax(T), realmin(Float32)
Out[5]:
In [6]:
T=BigFloat
eps(T), realmax(T), realmin(T), map(Int64,round(log10(1/eps(T))*log(10)/log(2)))
Out[6]:
We see that BigFloat has approximately 77 significant decimal digits (actually 256 bits) and very large exponents. This makes the format ideal for Greaffe's method.
Precision of BigFloat can be increased, but exponents do not change.
In [7]:
precision(BigFloat)
Out[7]:
In [8]:
setprecision(512)
eps(T), realmax(T)
Out[8]:
In [9]:
setprecision(256)
Out[9]:
Here is the function for the Graeffe's method. We also define small test polynomial with all real simple zeros.
In [10]:
using Polynomials
p=poly([1,2,3,4])
Out[10]:
In [11]:
roots(p)
Out[11]:
In [12]:
function Graeffe{T}(p::Poly{T},steps::Int64)
# map the polynomial to BigFloat
pbig=Poly(map(BigFloat,coeffs(p)))
px=Poly([zero(BigFloat),one(BigFloat)])
n=degree(p)
σ=map(BigFloat,2^steps)
for k=1:steps
peven=Poly(coeffs(pbig)[1:2:end])
podd=Poly(coeffs(pbig)[2:2:end])
pbig=peven^2-podd^2*px
end
# @show p[end]
y=Array{BigFloat}(n)
# Normalize if p is not monic
y[1]=-pbig[end-1]/pbig[end]
for k=2:n
y[k]=-pbig[end-k]/pbig[end-(k-1)]
end
# Extract the roots
for k=1:n
y[k]=exp(log(y[k])/σ)
end
# Return root in Float64
map(Float64,y)
end
Out[12]:
In [13]:
Graeffe(p,8)
Out[13]:
Now the Wilkinson's polynomial:
In [14]:
ω=poly(collect(one(BigFloat):20))
Out[14]:
In [15]:
Graeffe(ω,8)
Out[15]:
In [16]:
ans[2]
Out[16]:
In [17]:
Graeffe(ω,16)
Out[17]:
In [18]:
ans[2]
Out[18]:
We now generate the Chebyshev polynomial $T_{50}(x)$ using the three term recurence.
In [19]:
n=50
T0=Poly([BigInt(1)])
T1=Poly([0,1])
Tx=Poly([0,1])
for i=3:n+1
T=2*Tx*T1-T0
T0=T1
T1=T
end
In [20]:
T
Out[20]:
In [21]:
using Gadfly
In [22]:
f(x)=T(x)
Out[22]:
In [23]:
Gadfly.plot(f,-1,1)
Out[23]:
In order to use Graeffe's method, we need to shift $T$ to the right by one, so that all roots have simple moduli, that is, we compute $T(1-x)$:
In [24]:
Ts=T(Poly([BigFloat(1),-1]));
In [25]:
fs(x)=Ts(x)
Gadfly.plot(fs,0,2)
Out[25]:
In [26]:
# Computed roots, 16 steps are fine
y=Graeffe(Ts,16)-1
Out[26]:
In [27]:
# Exact roots
z=map(Float64,[cos((2*k-1)*pi/(2*n)) for k=1:n]);
In [28]:
# Relative error
maximum(abs,(z-y)./z)
Out[28]:
In [29]:
function mylu{T}(A1::Array{T}) # Strang, page 100
A=copy(A1)
n,m=size(A)
for k=1:n-1
for rho=k+1:n
A[rho,k]=A[rho,k]/A[k,k]
for l=k+1:n
A[rho,l]=A[rho,l]-A[rho,k]*A[k,l]
end
end
end
# We return L and U
L=tril(A,-1)
U=triu(A)
# This is the only difference for the block case
for i=1:maximum(size(L))
L[i,i]=one(L[1,1])
end
L,U
end
Out[29]:
In [30]:
A=rand(5,5)
Out[30]:
In [31]:
L,U=mylu(A)
Out[31]:
In [32]:
L*U-A
Out[32]:
In [33]:
L
Out[33]:
In [34]:
U
Out[34]:
We now try block-matrices. First, a small example:
In [35]:
# Try k,l=32,16 i k,l=64,8
k,l=2,4
Ab=[rand(k,k) for i=1:l, j=1:l]
Out[35]:
In [36]:
Ab[1,1]
Out[36]:
In [37]:
L,U=mylu(Ab)
Out[37]:
In [38]:
L[1,1]
Out[38]:
In [39]:
# Residual
R=L*U-Ab
Out[39]:
In [40]:
norm(R) # This is not defined
In [41]:
vecnorm(R)
Out[41]:
We need a convenience function to unblock the block-matrix:
In [42]:
unblock(A) = mapreduce(identity, hcat,
[mapreduce(identity, vcat, A[:,i]) for i = 1:size(A,2)])
Out[42]:
In [43]:
unblock(Ab)
Out[43]:
In [44]:
norm(unblock(L*U-Ab))
Out[44]:
We now compute timings and errors for bigger example:
In [45]:
# This is 512x512 matrix consisting of 16x16 blocks of dimension 32x32
k,l=32,16
Ab=[rand(k,k) for i=1:l, j=1:l]
# Unblocked version
A=unblock(Ab);
In [46]:
?lu
Out[46]:
In [48]:
# Built-in LAPACK function with pivoting
@time L,U,p=lu(A);
In [49]:
norm(L*U-A[p,:])
Out[49]:
In [50]:
# mylu() unblocked
@time L,U=mylu(A);
In [51]:
norm(L*U-A)
Out[51]:
In [53]:
# mylu() on a block-matrix - much faster, but NO pivoting
@time L,U=mylu(Ab);
In [54]:
norm(unblock(L*U-Ab))
Out[54]:
In [55]:
using Winston
In [56]:
k=20
n=20
E=Array{Any}(n,k)
# Unsymmetrix random uniform distribution
for i=1:k
A=rand(n,n)
E[:,i]=eigvals(A)
end
# We need this since plot cannot handle `Any`
E=map(eltype(E[1,1]),E)
Out[56]:
In [57]:
Winston.plot(E,"*")
Out[57]:
In [58]:
# Unsymmetric random normal distribution
E=Array{Any}(n,k)
for i=1:k
A=randn(n,n)
E[:,i]=eigvals(A)
end
# We need this for plot to work
E=map(eltype(E[1,1]),E)
Winston.plot(E,"*")
Out[58]:
In [59]:
# Symmetric random uniform distribution
E=Array{Any}(n,k)
for i=1:k
A=rand(n,n)
A=triu(A)+triu(A,1)'
E[:,i]=eigvals(A)
end
# We need this for plot to work
E=map(eltype(E[1,1]),E)
Winston.plot(E,"*")
Out[59]:
In [60]:
# Symmetric random normal distribution
E=Array{Any}(n,k)
for i=1:k
A=randn(n,n)
A=triu(A)+triu(A,1)'
E[:,i]=eigvals(A)
end
# We need this for plot to work
E=map(eltype(E[1,1]),E)
Winston.plot(E,"*")
Out[60]:
In [61]:
# Now the interactive partcxx
using Interact
In [62]:
@manipulate for k=10:30, n=10:30
E=Array{Any}(n,k)
for i=1:k
A=randn(n,n)
A=triu(A)+triu(A,1)'
E[:,i]=eigvals(A)
end
# We need this for plot to work
E=map(eltype(E[1,1]),E)
Winston.plot(E,"*")
end
Out[62]:
Mathematics is about spotting patterns! (Alan Edelman)
In [ ]: